Scintillation time dependence and pulse shape discrimination in liquid argon 
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Using a single-phase liquid argon detector with a signal yield of 4.85 photoelectrons per keV of 
electronic-equivalent recoil energy (keVee), we measure the scintillation time dependence of both 
electronic and nuclear recoils in liquid argon down to 5 keVee. We develop two methods of pulse 
shape discrimination to distinguish between electronic and nuclear recoils. Using one of these 
methods, we measure a background and statistics- limited level of electronic recoil contamination to 
be 7.6 X 10~^ between 52 and 110 keV of nuclear recoil energy (keVr) for a nuclear recoil acceptance 
of 50% with no nuclear recoil-like events above 62 keVr. Finally, we develop a maximum likelihood 
method of pulse shape discrimination based on the measured scintillation time dependence. 

PACS numbers: 61.25.Bi,29.40.Mc,95.35.-|-d 



I. INTRODUCTION 

Recent years have seen an increase in the number of ex- 
periments using noble Hquids as materials for detecting 
Weakly Interacting Massive Particles (WIMPs), a well 
motivated dark matter candidate The current best 
limit for the spin-independent WIMP-nucleon cross sec- 
tion for a 60-GeV WIMP mass is 4.6 x IQ-^^'cm^ set 
by the Cryogenic Dark Matter Search (CDMS) experi- 
ment 0]. Because many noble liquids have high scintil- 
lation yields, are easily purified of radioactive impurities, 
and are likely scalable to large masses with relative ease, 
they hold great promise for this application. The current 
best limit from a noble liquid detector is 8.8 x 10~^^ cm^ 
at 100 GeV, set by XENONIO 0]. Noble liquid detec- 
tors with larger target masses will likely improve on these 
limits. 

The key to the noble liquid dark matter detectors is 
discriminating between nuclear recoil events that consti- 
tute a WIMP signal and electronic recoil events that form 
the primary backgrounds. The XENON and ZEPLIN ex- 
periments are designed to collect both scintillation light 
and ionization from liquid xenon [J, Q . These are dual- 
phase detectors that use both scintillation light and ion- 
ization charge collection to discriminate between event 
classes, as nuclear recoils and electronic recoils produce 
different ratios of charge to light. Liquid argon is an 
attractive alternative to liquid xenon due to the lower 
cost of natural argon and its simpler purification require- 
ments. The WIMP Argon Programme (WARP) and 
Argon Dark Matter (ArDM) experiments employ dual- 
phase detectors that use liquid argon as the target [1, j?! . 

Alternatively, a single-phase detector collecting solely 
scintillation light might distinguish electronic and nuclear 



recoils using pulse shape discrimination (PSD). In 1977, 
Kubota et al. showed that the time dependence of scin- 
tillation light in liquid xenon and liquid argon is signifi- 
cantly different for heavy ionizers such as a particles and 
fission fragments when compared to light ionizers such 
as /3 decay and Compton-scattered electrons Q. This 
is because scintillation in liquid noble gases is produced 
by the decay of cxcimcrs that can exist in either sin- 
glet or triplet molecular states, which have very different 
lifetimes (Table HI . The slow scintillation light emitted 
by triplet molecules can be suppressed in intensity by 
destructive triplet interactions, primarily Penning ion- 
ization and electron-triplet spin exchange; it is believed 
that these reactions are stronger for high excitation den- 
sities such as those produced by nuclear recoils, causing 
the observed time dependences. Therefore, the relative 
amplitudes of the fast and slow components can be used 
to determine which type of excitation occurred for a given 
event. 

Singlet Lifetime (ns) Triplet Lifetime (ns) 



Ne 
Ar 
Xe 



< 18.2 ±0.2 
7.0 ± 1.0 
4.3 ±0.6 



14900 ± 300 
1600 ± 100 
22.0 ± 2.0 



'Electronic address: ' daniel.mckinsey@yale.edu| 



TABLE L Lifetimes of the sin glet and triplet states for neon, 
argon, and xenon excimers [gl. Il0||. 



Pulse shape discrimination based on the timing of scin- 
tillation light has been studied for use with several noble 
liquids. In liquid helium, PSD has been studied in or- 
der to separate electronic recoil events from ■^He(n,p)'^H 
events in the search for the permanent electric dipole 
moment of the neutron [ll|. PSD has also been used 
to suppress 7-ray backgrounds in liquid xenon [l^, [l^ . 
McKinsey and Coakley [13] pointed out that the much 
longer triplet lifetime in liquid neon should allow supe- 
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rior PSD, which has recently been verified experimen- 
tally Following this observation, Boulay and Hime 
recognized that the similar properties of liquid argon 
could in principle achieve PSD with part per billion lev- 
els of electronic recoil contamination (ERC) JJJ. ERC 
is defined to be the probability of incorrectly classifying 
an electronic recoil event as a nuclear recoil event given 
a particular level of nuclear recoil acceptance. An ERC 
of 10^^ or better is required to perform a competitive 
WIMP search using liquid argon due to the presence of 
the radioactive isotope ■^^Ar, which produces about 1 Bq 
per kg of atmospheric argon 

[M [O, El- The WARP 
collaboration has used scintillation timing in combina- 
tion with an ionization signal to reduce electronic re- 
coil backgrounds in liquid argon [To^. The Dark Mat- 
ter Experiment using Argon Pulse Shape Discrimination 
(DEAP) has demonstrated a background limited ERC 
of 5 X 10~^ using the DEAP-0 single phase detector for 
nuclear recoil energies above 1 MeV 20]. At a given en- 
ergy, ERC improves exponentially with scintillation light 
collection efficiency. For this reason, efficient scintillation 
light detection is the primary requirement for performing 
a sensitive WIMP search with negligible background at 
a suitably low energy threshold. 

In this paper we describe measurements of scintillation 
in liquid argon due to low-energy nuclear and electronic 
recoils in the energy range relevant to a WIMP dark mat- 
ter search. We measure the scintillation time dependence 
of liquid argon for both event classes. We develop two 
basic PSD methods and measure the level of discrimina- 
tion in our apparatus. Finally, we use the measured time 
dependence to develop a maximum likelihood method of 
PSD. 



II. EXPERIMENTAL DETAILS 
A. Detector Design 

The apparatus consists of a 3.14-liter active volume of 
liquid argon viewed by two 200-mm-diameter photomul- 
tiplier tubes (PMTs) [2l[, all contained within a stain- 
less steel vessel and vacuum Dewar. Figure [T] shows a 
schematic of the central volume and PMTs. The active 
region is defined by a Teflon cylinder 200 mm in diameter 
and 100 mm in height with two 3-mm-thick fused-silica 
windows enclosing the top and bottom. The PMTs are 
held in place by Teflon rings above and below the central 
volume and view the active region through the windows. 
They are powered by positive high voltage with a typical 
gain of approximately 4 x 10^. 

Because liquid argon scintillates in the ultraviolet 
(w 128 nm) ,22i]j slW inner surfaces of the Teflon and 
windows are coated with a thin film of tetraphenyl bu- 
tadiene (TPB) |23| by use of a Tectra Mini-Coater evap- 
oration system [24|[4J|. The TPB shifts the wavelength 
of the ultraviolet light to approximately 440 nm so that 
it may pass through the windows and be detected by 
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FIG. 1: (Color online) Schematic representation of the scin- 
tillation cell. 



the PMTs. Both windows are coated with (0.20 ± 0.01) 
mg/cm^ of TPB, while the Teflon cylinder is coated with 
(O.SOibO.Ol) mg/cm^. The Teflon cylinder, windows, and 
PMTs are all immersed directly in liquid argon and con- 
tained within a 25-cm-diameter by 91-cm-tall stainless 
steel vessel. 

The stainless steel vessel is in turn housed inside a vac- 
uum Dewar, and argon gas is introduced into the system 
though a tube on the top of the Dewar. The argon is liq- 
uefled in a copper cell mounted to the end of a pulse-tube 
refrigerator |25j inside the Dewar before flowing through 
a tube to the stainless steel vessel. All components that 
come into contact with the gas or liquid are baked to at 
least 60°C, and the ultra-high-purity argon gas (99^999%) 
is passed through a heated gas-purification getter [26] be- 
fore entering the vessel. In addition, the argon is contin- 
ually circulated through the getter and reliquefied at a 
rate of at least 2.0 standard liters per minute (slpm) to 
ensure that high purity is maintained. The stability of 
the system is discussed further in Sec. Ill CI 

The data acquisition system is custom-built around 
VME-bus waveform digitizers (WFDs); a sample WFD 
trace from a scintillation event in argon can be seen in 
Fig. [21 The PMT signals from the detector are divided 
three ways by a linear fan out with two copies of the 
signal sent to the WFDs and one sent to a triggering 
system. Each WFD has four channels that record eight- 
bit samples at 500 MHz. These samples are stored in 
a separate programmable-length memory buffer for each 
channel. For all data presented here, the record length 
is set to 26 /is. The two copies of the PMT waveforms 
are recorded separately at unity gain and at an atten- 
uation of ten to increase the effective dynamic range of 
the eight-bit digitization. The buffer is continually filled 
but saved to disk only when the triggering system regis- 
ters a fraction of a photoelectron in both PMTs within 
a 100-ns coincidence window. Once a trigger has been 
registered, the DAQ records the event for 22 fj,s, leaving 
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an additional 4 fj,s of baseline presamples in the data (see 
Fig. [2) . The data are read to a computer via fiber-optic 
cable, and after the computer has recorded all 26 /iS of 
data, it resets the system for the next event. The data 
collected by the DAQ software are saved in a ROOT- 
based file structure 12711. 




20 25 
Time (|Xs) 

FIG. 2: Example of an electronic recoil event from a single 
PMT, digitized by the 8-bit WFD, sampling at 500 MHz. 
Four microseconds of presamples are recorded to measure the 
baseline. 



B. Data Collection 

All data collected are processed in software. First, for 
each PMT, the two gain scales are combined into a single 
waveform. Three /zs of baseline presamples are averaged 
to obtain a baseline and baseline root-mean-square. The 
baseline is then subtracted from the trace. Ideally, we 
would count single photoelectrons in a pulse, but since 
we detect many photoelectrons that produce signals over- 
lapping in time, we integrate the trace in order to deter- 
mine the total number of photoelectrons. To mitigate 
the integration noise, we restrict the range of the inte- 
gral to 50-ns regions in which the trace voltage crosses 
a threshold of 4 times the baseline root-mean-square (a 
threshold of approximately 2/5 the height of a photoelec- 
tron). This method is a hybrid of single photoelectron 
counting and pulse integration. 

We then apply three cuts to all data. One cut removes 
events in which either PMT saturates due to excessive 
light exposure by rejecting events above an experimen- 
tally determined threshold of approximately 2000 times 
the single photoelectron pulse area. A second cut re- 
moves events for which the trigger time (defined as the 
time at which the voltage rises above 20% of its maxi- 
mum value) differs by greater than 20 ns between the two 
PMTs. The third cut is designed to eliminate events that 
produce light in the windows or the glass of the PMTs. 
An asymmetry parameter A is defined as 



where St and Sb are the signal areas in the top and 
bottom PMTs. For most data, we require —0.275 < A < 
0.375, due to a slightly larger gain in the top PMT. This 
cut is relaxed to \ A\ < 0.4 to improve the statistics of the 
nuclear recoil data described below. 

We use a 10-/xCi ^^Na source to produce electronic 
recoils. In 90% of ^^Na decays, a positron is emitted 
that immediately annihilates in the surrounding materi- 
als to produce 511-keV 7 rays with equal and opposite 
momenta, or "back-to-back". We use the second 7 ray 
to tag electronic recoil events in the liquid by triggering 
on a coincidence within a 100-ns window between the 
PMTs in the argon and a Nal crystal scintillator placed 
back-to-back with our apparatus. This event tagging re- 
duces backgrounds in our data from other radioactive 
decays and cosmic rays. To further decrease neutron 
backgrounds, we place one layer of water-filled containers 
above and around the sides of the dewar. These contain- 
ers are cubes of side 30 cm in length and hold 20 liters 
of water. 

In addition to the three universal cuts, we apply two 
additional cuts to ensure data purity. In software we 
narrow the coincidence window between the liquid argon 
PMTs and the Nal crystal to 30 ns. We also make an 
energy cut to the NaT energy spectrum to select 511- 
keV events in the Nal crystal. In the liquid argon, the 
majority of the 511-keV 7s Compton scatter, producing 
events with a continuum of deposited energies. 

To investigate the detector response to nuclear recoils, 
we use a portable deuterium-deuterium neutron gener- 
ator [2^ as a source of 2.8-MeV neutrons and a PMT 
viewing BC501A organic scintillator as a secondary de- 
tector. Both the generator and the organic scintillator 
are placed approximately 1.63 m from the center of the 
active volume. We require a detection of a scintillation 
event in the liquid argon, followed within 200 ns by an 
event in the organic scintillator. The experimental setup 
can be seen schematically in Fig. [3] By changing the 
scattering angle 9, we can choose the energy of the nu- 
clear recoils observed in the liquid, Erec, using simple 
kinematics; 



Ejre.c- 



(T+mF 



[H-M-cos2(6i) 



3(61) ^Af2 -f cos2(6') - 1] , 



(2) 



A 



St — Sb 
St + Sb' 



(1) 



where Ein is the incident neutron energy and M is the 
atomic mass of the target. 

Since we know the energy of the neutron, we can cal- 
culate the time-of-flight of a neutron that scatters in the 
liquid argon and the organic scintillator. We then ap- 
ply a time-of-flight cut between the liquid argon cell and 
the organic scintillator to distinguish between 7 rays and 
neutrons, as the neutron generator produces both. In 
general, the time-of-flight cut requires an event to occur 
in the organic scintillator 60-90 ns after the event in the 
liquid argon. In addition, we make a PSD cut in the 
organic scintillator data to further eliminate electronic 
recoils. 
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FIG. 3: (Color online) Schematic of the neutron scattering 
setup. 



measure the signal yield of the detector to be 4.85 pho- 
toelectrons per keV electron equivalent (keVee), where 
keVee refers to the amount of energy deposited by an 
electronic recoil. The response of the detector was stable 
to within 5% during the four months of data acquisition. 




100 200 300 400 500 600 700 800 
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We also collect background data by looking at events 
in the liquid argon with no external source present. This 
background data provides an estimate of the accidental 
background rate that may be contaminating the 7-ray 
data sets, and it will be discussed further in Sec. IIII Bl 



C. Detector Calibration 

We use a 10-/iCi sealed ^^Co source for daily calibra- 
tions. This source produces 122-, 137-, and 14.4-keV 7 
rays, with branching ratios of 86%, 11% and 9%, respec- 
tively. Any scintillation event in liquid argon produces 
a significant triplet component; since this component is 
spread out over many fis, it appears in the signal as many 
single photoelectrons well separated in time (for exam- 
ple, there are a number of single photoelectrons that ap- 
pear after 7 fj,s in Fig. Therefore, we measure the 
gain of the PMTs using the ^^Co source by selecting sin- 
gle photoelectrons from the tail end of each pulse. The 
PMT traces are divided into 75-ns regions centered on 
times at which the trace crossed an experimentally de- 
termined threshold of roughly 1/3 of a photoelectron. 
These regions are then integrated to obtain the single 
photoelectron pulse area. The typical gain for the PMTs 
is approximately 4 x 10^. 

Figure m shows an example ^''Co spectrum along with 
a simulation done with the Reactor Analysis Tool, a 
toolkit of Geant4 developed by the Braidwood collab- 
oration [2^, [sO]. From the simulation, we find that the 
position of the primary peak is dominated by the 122- 
keV 7 photoabsorption process. The simulation parame- 
ters describing absorption and reflection of the materials 
in the detector are tuned to match the observed signal 
yield, and the spectral shapes line up nicely. 

By comparing the integrated signal corresponding to 
the 122-keV peak to that of a single photoelectron, we 



FIG. 4: (Color online) Example spectrum of ^^Co data, along 
with a simulation done using the Reactor Analysis Tool [3(^ . 

We use the 122-keV ^^Co peak to provide a daily en- 
ergy calibration. To check the quality of that calibration, 
we use the 511-keV 7 rays produced by the ^^Na source 
as a second point of reference. When calibrated using the 
^^Co source and the assumption that the signal scales lin- 
early with deposited energy, the 511-keV absorption line 
appears in the ^^Na spectrum as expected to within 1%. 
In addition, simulations of 511-keV 7 rays are consistent 
with the data. 

Impurities in the detector can build up over time via 
outgassing. These impurities can quench argon excimers 
or absorb emitted UV photons, which would lead to a de- 
crease in light yield. Additionally, one would observe a 
decrease in the triplet molecule lifetime. Work by Himi et 
al. suggests that an impurity level of 0.5 atoms of ni- 
trogen per 10^ atoms of argon in the liquid could de- 
crease the observed triplet lifetime by as much as 0.1 
/iS 31]. Further experiments quantifying the reduction of 
the triplet lifetime due to nitrogen and oxygen impurites 
have recently been performed by the WARP collabora- 
tion [33, [3^ • To avoid signal degradation, we continually 
circulate the argon through a getter before reliqucfying 
back into the detector. We use daily measurements of 
both the light yield and the triplet lifetime to monitor 
the purity level. 

We measure the light yield in the manner described 
earlier in this section, and we use the same ^^Co data to 
measure the triplet lifetime. First, we select events in the 
122-keV peak to make sure we use a similar data set for 
each individual measurement. The top and bottom PMT 
traces are normalized by the size of the single photoelec- 
tron and summed together. We align each pulse based on 
its estimated trigger time, defined as the time at which 
the trace first crosses 20% of its maximum value. At this 
trigger time, the relative time for each pulse is t = 0. 
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Between 5000 and 10000 traces are averaged and the fol- 
lowing model is fit to the average trace between 1 and 7 
^s from the trigger: 



<V{t) >^ Aexp{-t/Ti) + B, (3) 



where < V{t) > is the expected trace, A is a normaliza- 
tion factor, Ti is the triplet lifetime and B is an additional 
baseline term that helps stabilize the fit over a range of 
fit windows. An example fit is shown in Fig. [5l 




012345678 
Time from trigger (jls) 



FIG. 5: (Color online) An example fit of a single exponential 
to an average trace to measure the long time constant in liquid 
argon. 



We find two significant sources of systematic error in 
the fit parameters stemming from the voltage applied to 
the PMTs and the choice of fit window. Although we 
do not fully understand the PMT voltage effect, we esti- 
mate the systematic uncertainty to be 33 ns by changing 
the PMT voltages by ±75 V (equivalent to dividing or 
multiplying the gain by 2). We do not fit the data before 
1 fis or after 9 /zs, as the fit becomes less reliable due 
to contamination from the fast component and baseline 
noise, respectively. We estimate the systematic uncer- 
tainty associated with the fit window to be about 35 ns 
by both varying the end time by ±2 fis and by fitting the 
data within 5 ^s windows ending at 6, 7 and 8 /is from 
the trigger. We combine the two sources of error into a 
single systematic uncertainty estimate of 50 ns. 

During the four months of operation, we mea- 
sure the long time constant and signal yield to be 
(1463 ± 5stat ± 50sys) ns and (4.85 ± 0.08) photo- 
electrons/keVee. The uncertainty on the signal yield is 
statistical, and varying the PMT voltage by ±75 V has 
no apparent effect. For a fixed fit window and PMT volt- 
age, the long time constant is stable to within 1%, while 
the signal yield is stable to within 5%. 



III. EXPERIMENTAL RESULTS 

A. Detection Time PDF Model 

We measure the time dependence of scintillation light 
produced by electronic and nuclear recoil scattering 
events in our detector. For an event with energy de- 
posited at time to, we model the temporal probability 
density functions (PDFs) for the emission times of scin- 
tillation photons as the weighted sum or mixture of two 
exponential PDFs: 

fit -to) =q g{t-to,Ti) + {l-q) g{t-to,Ts), (4) 

where 




Ti and Ts are the long and short time constants for both 
nuclear and electronic recoil event classes at any given 
energy, and the probability parameter q takes different 
values for the two classes of events. 

At each PMT, detected scintillation photons yield pho- 
toelectrons that produce observed voltage traces. For 
the small detector in this experiment, scintillation tran- 
sit times are negligible. We assume that the duration of 
energy deposition and excimer formation is instantaneous 
compared to the time scales relevant to scintillation light 
emission. Hence, we model the expected voltage trace 
as a convolution of the impulse response function of the 
PMT and the PDF model for the emission of scintillation 
photons: 

/>oo 

<V{t)>(x hv{t-s)f{s-t^)ds, (6) 

Js=t, 

where hy (t) is the impulse response function of the PMT 
and is an additional model parameter that relates the 
energy deposit time to to the relative time scale we as- 
sociate with our measurement of hv{t). We estimate the 
impulse response function by averaging single photoelec- 
tron events observed in our calibration data. Neglecting 
additive noise and other instrumental systematic errors, 
the integral of the voltage trace is proportional to the 
number of photoelectrons produced by the event. 

The data are divided into 15 non-uniform bins by pho- 
toelectron number, with the smallest signal bin includ- 
ing events consisting of 20-24 photoelectrons and the 
largest bin including events consisting of 240-279 pho- 
toelectrons. These bins define the region of interest. 
For each photoelectron bin, we generate template traces 
for both electronic and nuclear recoils by averaging all 
"tagged" events of a given type. The traces are aver- 
aged in the same manner as described in Sec. Ill CI and 
normalized. 

We use these average voltage traces to determine the 
model parameters: t;, Ts, q„uciear and qeiectronic- The ob- 
served and predicted fraction of a normalized trace in the 



6 



ith time bin are called Pm{i) and p{i)^ respectively. We 
obtain the model parameters by minimizing the squared 
Matusita distance [s^, [H, [IBl between pm and p: 

\P-Pvi\m = ^WPrn{i) - \f¥(^)f, (7) 

i 

where negative values oipm are set to 0. For each bin, we 
fit our model to normalized mean voltage trace data in 
a time window that ends about 6800 ns after the trigger 
time. We determine to be 30 ns before the trigger time, 
chosen to minimize the value of Eq.[7l Table Ull shows the 
estimated model parameters for each photoelectron bin. 
Figure [6] shows the model prediction along with data for 
the 80-99 photoelectron bin. 

Between 80 and 300 ns, there is a feature in both event 
classes that is not well predicted by the model. A scintil- 
lation component that decays as approximately has 
been observed in liquid helium, attributed to diffusion- 
dominated excimer-excimer destruction flU]. When we 
fit a model including the full form of this component as 
described in the parameter estimates are not well 
determined. Previous observations of scintillation in ar- 
gon have observed an intermediate exponential compo- 
nent with a decay time of 20-40 ns @ . Unfortunately, a 
model that includes a third exponential component nei- 
ther returns stable model parameters as a function of en- 
ergy nor accurately predicts the trace behavior between 
80 and 300 ns, so we prefer the two-component model. 
We assume that the sharp bump localized at 150 ns is 
caused by the cabling and electronics. 



Bin (pe) 


Ti (ns) 


Ts (ns) 


Qjiuclear 


^electronic 


20- 


-24 


1634 ± 150 


9±3 


0.378 ±0.011 


0.523 ±0.015 


25- 


-29 


1535 ± 128 


10 ±3 


0.382 ± 0.011 


0.573 ±0.015 


30- 


-34 


1478 ± 107 


10 ±3 


0.357 ±0.010 


0.601 ±0.014 


35- 


-39 


1455 ± 102 


11±3 


0.353 ±0.010 


0.627 ±0.014 


40- 


-49 


1461 ± 96 


12 ± 3 


0.344 ±0.010 


0.658 ±0.014 


50- 


-59 


1459 ± 92 


12 ± 3 


0.327 ± 0.009 


0.681 ±0.015 


60- 


-69 


1439 ± 89 


12 ± 3 


0.315 ±0.010 


0.699 ±0.015 


70- 


-79 


1448 ± 89 


13 ±3 


0.309 ±0.010 


0.710 ±0.015 


80- 


-99 


1447 ± 85 


13 ±3 


0.298 ±0.010 


0.721 ±0.015 


100- 


-119 


1452 ± 84 


13 ±3 


0.289 ±0.010 


0.733 ±0.015 


120- 


-139 


1447 ± 84 


13 ±3 


0.284 ±0.011 


0.741 ±0.016 


140- 


-159 


1446 ± 84 


14 ±3 


0.278 ±0.012 


0.747 ± 0.016 


160- 


-199 


1450 ± 84 


14 ±3 


0.272 ±0.013 


0.752 ± 0.016 


200- 


-239 


1460 ± 84 


15 ±3 


0.265 ±0.015 


0.760 ± 0.016 


240- 


-279 


1467 ± 84 


15 ±3 


0.258 ±0.018 


0.764 ± 0.016 



TABLE 11: Estimated model parameters and 1-sigma uncer- 
tainties for each photoelectron bin. The systematic errors 
described in the text are the dominant source of error. 



We determine 1-sigma random uncertainties for the 
model parameters with a nonparametric bootstrap re- 
sampling scheme [131 • A rigorous quantification of sys- 
tematic uncertainty arising from the fit window, the 




Relative time (ns) 




Relative time (ns) 



FIG. 6: (Color online) Observed and predicted mean volt- 
age traces for nuclear and electronic recoil events of 80 to 99 
photoelectrons. 



voltage applied to the PMT, and variation in the re- 
sponse of each individual PMT is difficult. Nonetheless, 
we get approximate estimates of systematic uncertain- 
ties from these sources, and these systematics are the 
dominant source of error. We estimate 1-sigma system- 
atic uncertainties in the parameters Ti (for i — I, s) as 

0.5 X {n ,max — Ti^rnin)^ where Ti^rnax and Ti^min BIG the 

maximum and minimum values taken by Ti for different 
choices of time windows. We refit the model to data using 
time windows of approximately 5000 ns and 8600 ns from 
the trigger time, and vary the offset parameter by ± 
2 ns. We also check for differences between the response 
of each PMT in five different photoelectron bins. We in- 
clude the voltage-dependent errors described in Sec. Ill CI 
in the systematic uncertainty of ti , and we estimate the 
voltage-dependent systematic variation of the g-values. 
We find that the time constants are most influenced by 
choice of fit window, while the g-values are most infiu- 
enced by PMT effects. We combine all sources of error 
to obtain the estimates of systematic uncertainty shown 
in Table ni 

As a consistency check on the quality of fit, for each 
photoelectron bin we use our model to predict the prompt 
fraction, a discrimination statistic described in detail in 
the next section. We compute the difference between 
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the predicted and observed prompt fractions, and the 
root-mean-square values of this prediction error across 
all bins are respectively 0.007 and 0.003 for the nuclear 
and electronic recoil event classes. The fractional root- 
mean-square value of this prediction error is about 1% 
for both event classes. 

B. Prompt Fraction Method 

The prompt fraction method is a simple approach to 
pulse shape discrimination. For each trace, we define the 
prompt fraction fp as 

_ 4 vm 

where V{t) is the voltage trace from the PMT, ^ is 
an integration time determined to optimize the ERC, 
Ti — to — 50 ns, Tf — to + 9 ^s, and to is the trigger 
time as defined in Sec. Ill Bl The measured discrimina- 
tion does not significantly improve by extending Tf to 
20 /xs. Fig. [7] shows a scatter plot of fp versus energy 
for both electronic and nuclear recoils. The two popula- 
tions of events represent the tagged data remaining in the 
neutron generator data set and the ^^Na data set after 
the selection cuts described in Sec. Ill Bi have been made. 
We choose ^ = 90 ns by estimating the ERC based on a 
simple Gaussian model for values of ^ from 50 to 250 ns 
over a variety of different photoelectron bins. Although 
the value of ^ has only a weak effect on the predicted 
discrimination, a choice of ^ = 90 ns provides the best 
results across the widest range of energies in the region of 
interest, and that value is used for the analysis presented 
in this paper. 




F .... I .... I .... I .... I I ... I 
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Energy (keVee) 

FIG. 7: (Color online) A scatter plot of fp vs. energy for 
tagged electronic and nuclear recoils, where ^ = 90 ns. 

For energy bins of width 1 keVee between 5 and 32 
keVee, we form histograms of the electronic and nuclear 
recoil fp statistics. To estimate the expected value of fp, 
fp, we fit a Gaussian function to the empirical distribu- 
tions. In Table IIIII and Fig. [51 we present the estimated 



mean fp for both classes of events in the energy range 
of interest. We estimate the systematic uncertainty on 
the values in Table IIIII to be 3%. This uncertainty esti- 
mate comes from changing the PMT voltage by — 75 V 
and from variations in the measured signals between the 
two PMTs. The mean values of the fp distributions for 
the two event classes are closer at low energies than at 
high energies, possibly because dE/dx for nuclear recoils 
decreases at low energies while increasing for electronic 
recoils. Therefore, the PSD improves at higher energies 
both because of increased photoelectron statistics and be- 
cause of increased separation between the mean fp val- 
ues. 



En6rgy (k6"V66) fp, electronic fp, nuclear 

5- 6 0.391 ± 0.012 0.566 ± 0.018 

6- 7 0.376 ± 0.011 0.595 ± 0.018 

7- 8 0.361 ± 0.011 0.607 ± 0.019 

8- 9 0.349 ± 0.011 0.625 ± 0.019 

9- 10 0.339 ± 0.010 0.638 ± 0.020 

10- 11 0.334 ± 0.010 0.640 ± 0.020 

11- 12 0.328 ± 0.010 0.649 ± 0.020 

12- 13 0.322 ± 0.010 0.663 ± 0.020 

13- 14 0.319 ± 0.010 0.658 ±0.020 

14- 15 0.314 ± 0.009 0.675 ± 0.020 

15- 16 0.311 ± 0.009 0.683 ± 0.021 

16- 17 0.309 ± 0.009 0.678 ± 0.021 

17- 18 0.304 ± 0.009 0.685 ± 0.021 

18- 19 0.302 ± 0.009 0.682 ± 0.021 

19- 20 0.299 ± 0.009 0.684 ± 0.021 

20- 21 0.297 ± 0.009 0.690 ± 0.021 

21- 22 0.295 ± 0.009 0.695 ± 0.021 

22- 23 0.292 ± 0.009 0.699 ± 0.021 

23- 24 0.290 ± 0.009 0.690 ± 0.021 

24- 25 0.288 ± 0.009 0.688 ± 0.021 

25- 26 0.289 ± 0.009 0.695 ± 0.021 

26- 27 0.288 ± 0.009 0.696 ± 0.021 

27- 28 0.285 ± 0.009 0.696 ± 0.021 

28- 29 0.284 ± 0.009 0.701 ± 0.021 

29- 30 0.283 ± 0.009 0.708 ± 0.021 

30- 31 0.281 ± 0.009 0.701 ± 0.021 

31- 32 0.282 ± 0.009 0.689 ± 0.021 



TABLE III: This table presents estimated mean values of fp 
versus energy, where ^ = 90 ns. The main sources of un- 
certainty are systematic, stemming from voltage effects and 
differences in the measured signals between the two PMTs. 

We estimate the ERC as the number of tagged elec- 
tronic recoil events with fp > /p, nuclear divided by the 
total number of electronic recoil events, where /p, nuclear is 
the estimated mean fp for nuclear recoils of that energy. 
This restriction sets a nuclear recoil acceptance level of 
approximately 50%. Since the shielding, coincidence, and 
timing cuts do not eliminate all neutron backgrounds in 
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FIG. 8: (Color online) The estimated mean fp versus energy 
for both event classes where ^ = 90 ns. 



the detector, a background estimation iVbg is made by 
measuring the rate of background neutrons i?bgn and as- 
suming that the background is dominated by these neu- 
trons hitting in accidental coincidence with 7 rays (Rj) 
in the liquid scintillator during the time allowed by the 
time-of- flight cut: 



i?b^n X R^X TOF X Ta, 



where Ta is the acquisition time of the data. To enable 
comparison with the measured ERC, we divide Nbg for 
each energy bin by the total number of electronic recoil 
events in that bin. 

Figure [9] shows the ERC observed using the prompt 
fraction method. We also plot the background estima- 
tion, two PSD projections based on the statistical model 
described below, and the ERC observed by applying a 
multibin method of PSD described in the next section. 
We convert the energy axis in Fig. [9] to keV of nuclear 
recoil energy (keVr) from keVee by dividing all electron 
equivalent energies by a constant nuclear recoil scintilla- 
tion efficiency of 0.29. This value was obtained by mea- 
surements using the same apparatus described in this pa- 
per and will be discussed in an upcoming publication [38| . 
We present the PSD results in keVr because that is the 
unit of interest for a dark matter detector. Using the 
prompt fraction method, for a nuclear recoil acceptance 
level of approximately 50%, we measure a background- 
and statistics-limited level of ERC in our detector of 
8.5 X 10~^ between 52 and 110 keVr (11 contamination 
events). We observe no nuclear recoil- like events above 
69 keVr. For comparison, there is an uncorrelated neu- 
tron background rate of 6 mHz between 69 and 110 
keVr, corresponding to 0.25 expected background counts 
in that energy range. 

Follo wing work done by members of the DEAP collab- 
oration |39l . l40l | , we model our estimates of the number 
of photoelectrons in the prompt and late time windows, 
Np and Ni, as normally distributed, independent ran- 



FIG. 9: (Color online) Measured electronic recoil contami- 
nation obtained by using the prompt fraction method and 
the multi-bin method. Also shown: the background estima- 
tion and model predictions described in the text. We include 
a model with the additional noise set to for comparison. 
There were no contamination events above 69 keVr observed 
by use of either method. The energy axis has been scaled from 
keVee to keVr by use of a constant nuclear recoil scintillation 
efficiency of 0.29 as discussed in the text. 



dom variables with means fip and and variances cTp 
and af. The estimated total number of photoelectrons, 
— Nj) + Ni, is also a random variable, with mean 



+ Hi and variance (7^^^ =07, + af. We express 



Ntot 
Mtot 

fip and Hi in terms of /Ztot and fp and we decompose the 
variances into two components: 



2 

Op ^ 



/pMtot 

(1 - /p)A*tot 



2 



(10) 



where cTp. add and i7;.add represent additional sources of 
random variability beyond what we expect from Poisson 
counting statistics (for example, integration noise). 

Hinkley [4lj has described in detail the probability den- 
sity function of the ratio of two normally distributed, cor- 
related random variables. For simplicity, we present here 
an approximation [given by Eq. 9 of Ref. [4l|] to the PDF 
of /p = iVp/7Vtot: 



IT, fj,j,x+a^tii (l-x) 
V2;r(o-fa;2+cr2(l-a;)2)3/2 ^ 

(fiix-fj.p(l-x)^) 1 
2{afx2+al{l-xy)i 



exp[ 



(11) 



where we have used the fact that the correlation, p, 
tween Np and A'tot is 



be- 



(12) 
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For the analysis below, we use the exact PDF [Eq. 1 
in Ref. (41l|]. although in practice the approximation is 
extremely good down to the lowest energy bin examined. 

We fit the electronic recoil data with the statistical 
model in each energy bin by fixing TVtot according to our 
measured light yield and treating fp, cTp^add, and (T;,add as 
free parameters. We assume that the statistical distribu- 
tion of fp does not strongly depend on the value of iVtot 
for events of the same energy. For any particular energy 
deposit, A^tot is a random variable, and events due to 
many different energy deposits can contribute to prompt 
ratio data in any one bin in A'tot space. The probabil- 
ity density function for prompt ratio data is a mixture 
of energy dependent PDFs, and we neglect this energy 
blurring effect. Simple Monte Carlo studies suggest that 
the ratio-of-Gaussians model breaks down for idealized 
Gaussian data due to the constraint on TVtot resulting 
from the binning of data, and this effect has not been 
taken into account in our analysis. 

We use the statistical model to estimate the expected 
fraction of electronic recoils that are misclassified as nu- 
clear recoils: 

ERC= f 9f^{x)dx. (13) 

Here, we choose r] = /p, nuclear to set the nuclear recoil 
acceptance level to approximately 50%, and we choose 
the parameters of gf^ according to the fits to the elec- 
tronic recoil fp distribution. Fig. [5] shows the predicted 
ERG versus energy according to this model. We also plot 
the idealized case where CT;,add and (ip.add are set to 0. In 
general, we expect that the normal distribution model for 
Np and Ni is an approximation for any energy of interest. 
For the idealized case where fip^add = ci^add = 0, Np and 
Ni would be Poisson random variables rather than Gaus- 
sian random variables. Hence, the ERG predicted by the 
ratio-of-Gaussians model for the idealized case should be 
interpreted with caution particularly at lower energies 
where the accuracy of a normal distribution model for a 
Poisson random variable can be very poor. 

Figure [TUl shows an example of the model fit for 14-15 
keVee and 30-31 keVee electronic recoil events. There 
is a deviation from the model at low fp values that we 
attribute to pile-up and noise triggers. There is also an 
excess of events in the high fp region. This excess might 
be caused by edge effects in our detector coming from 7 
tracks going into the walls or the TPB layer, producing 
extra prompt light. A larger detector with position re- 
construction capability might be able to eliminate such 
edge effects. There could be some unknown phenomena 
at work in the production and decay of argon molecular 
states, yielding a small fraction of events with anoma- 
lously large fp values. A third possibility is approxima- 
tion error stemming from incorrect model assumptions, 
such as the Gaussianity of Np and Ni or the effect of data 
binning. These possible effects can be better investigated 
by detectors with improved neutron shielding, and fu- 



ture studies will require more work to better calibrate the 
ERG predictions made by the ratio-of-Gaussians model. 
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FIG. 10: (Color online) Projections of the electronic recoil 
data from Fig. [7] onto the y-axis for 14-15 keVee (top) and 
30-31 keVee events (bottom), fitted by Eq. [TT]as discussed in 
the text. 

Figure [H] shows the fitted values of the variances, 
af = m + erf and = fip + a^^^^^ versus the es- 
timated number of photoelectrons in /i/ and /ip. The 
ideal case of zero additional noise (cr^ = /i) has been 
shown for comparison. The relationship between the es- 
timated variance parameters and the corresponding esti- 
mates of the number of photoelectrons is approximately 
linear, suggesting that cr^^add and are proportional 

to /i; and /ip. Straight line fits to the late and prompt 
variances have slopes of 2.2 and 1.3, respectively. Due to 
the various uncertainties in the fit such as that associ- 
ated with data binning, these numbers do not represent 
a rigorous estimate of the amount of noise in the late and 
prompt distributions; however, they might be used as a 
point of comparison with future detectors in trying to re- 
duce the overall noise in the system. Possible sources of 
this additional noise include integration noise, the widths 
of the single photoelectron spectra, and variability in the 
prompt window size arising from uncertainty in the de- 
termined trigger position. To improve PSD, future ex- 
periments will need to reduce the size of this additional 
noise to approach the ideal case. 
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FIG. 11: (Color online) Estimated af and (Tp parameters from 
the statistical model plotted against the estimated number of 
photoelectrons in the late and prompt components, respec- 
tively. Each distribution is fitted by a straight line, and the 
ideal case where — fi is also shown. 



C. Multibin Method 

The prompt fraction method of PSD is based on bin- 
ning the voltage trace into two time bins. We generalize 
this approach by representing a normalized voltage trace 
as a X L dimensional matrix. In this representation, we 
categorize the data by the number of photoelectrons in 
the event, and L refers to this division of data into photo- 
electron bins. We also partition the voltage trace into K 
time bins. We choose K = 10 for ease of computation, 
but there may be a better choice of K. With the ex- 
ception of the four smallest signal bins, which have been 
combined to form bins of 20-29 and 30-39 photoelec- 
trons, we use the same average voltage traces obtained 
in Sec. IIII Al as templates. 

We partition each template into K — 10 time bins. 
The upper and lower endpoints of each time bin are se- 
lected so that the fractions of the template trace for 80- 
99 photoelectron electronic recoil events that falls in each 
time bin are approximately equal. The initial bin starts 
50 ns before the trigger, and the endpoints of each bin 
are as follows, measured in nanoseconds from the trig- 
ger: 8, 18, 56, 200, 440, 750, 1180, 1800, 2950, and 8000. 
We do not adjust these endpoints for different photoelec- 
tron number but use the 80-99 photoelectron based time 
binning scheme for all cases. 

The kth component of the normalized template for an 
event that falls in the Ith photoelectron bin for the nu- 
clear and electronic recoil classes is denoted p„(fc,Z) or 
Pe{k, I), respectively. For instance, Pn{l, 1) represents the 
fraction of the 20-29 photoelectron nuclear recoil tem- 
plate trace that falls in the first time bin (between 50 ns 
before the trigger and 8 ns after the trigger). 

To assign an event to either the nuclear recoil or elec- 
tronic recoil class, we first compute a discrimination 
statistic for the event. We motivate a discrimination 



statistic based on analysis of an idealized experiment in 
which we could observe the absolute detection time of 
each photoelectron without error. In this ideal case, for 
a fixed number of detected photoelectrons, the observed 
number of photoelectrons in the time bins is a multino- 
mial random variable. 

Given that the fraction of detected photoelectrons in 
the fcth bin is pm{k,l), the multinomial log-likelihood 
statistics for the nuclear and electronic recoil classes are 

K L 

6w pmik, I) In pn{k, I') + const (14) 

fc=i 1=1 

and 

K L 

Sii'Pm{k, I) \npe{k, I') + const, (15) 

k=l 1=1 

respectively Here, A'tot is the total number of de- 

tected photoelectrons, and Pn{k,l) and Pe{k,l) are the 
expected values of the fraction of detected photoelectrons 
in the fcth time bin for the nuclear and electronic recoil 
classes, respectively. For this idealized case, a natural 
choice for the discrimination statistic is the log-likelihood 
ratio statistic, InRm'- 

lni?„ = lnr, -lnr„. (16) 

In our experiment, the observed data is not multino- 
mial, primarily because we observe a noisy voltage wave- 
form rather than discrete detection times. Nonetheless, 
we compute a discrimination statistic using Eqns. fMHTBl 
where we estimate iVtot and pmik,l) from any voltage 
trace of interest, and determine pe{k,l) and Pn{k,l) as 
described earlier. 

Figure [T^] shows a scatter plot of In Rm versus energy, 
analogous to Fig. [71 From this point the analysis parallels 
the prompt fraction method, as we form histograms of 
InRrn by energy bin and fit a Gaussian function to the 
observed In statistics to estimate the mean of lni?,„. 
Fig. [13] shows an example of the fitted projections for 
14-15 keVee and 30-31 keVee events. 

For each energy bin, we estimate the mean values of 
In R„i for nuclear recoils based on the Gaussian fits to the 
observed distributions, and this estimated mean value de- 
termines an approximate 50% nuclear acceptance thresh- 
old. We then determine the fraction of events in the 
tagged electronic recoil data set that have discrimination 
statistics less than this mean value to determine the ob- 
served level of ERG using the multi-bin method. If we 
determine the 50% nuclear acceptance threshold by find- 
ing the median values of the nuclear recoil distributions, 
the observed ERG is not significantly affected. Figure [5) 
in Sec. IIII B] shows the ERG using the multibin method. 
The mtiltibin method outperforms the prompt fraction 
method by as much as an order of magnitude. For a nu- 
clear recoil acceptance of approximately 50%, we have 
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FIG. 12: (Color online) A scatter plot of In Rm vs. energy for 
both electronic and nuclear recoils. 




Log-likelihood statistic, In S, 
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FIG. 13: (Color online) Projections of Fig. ll2lonto the y-axis 
for 14-15 keVee (top) and 30-31 keVee (bottom) events, with 
Gaussian fits to both the electronic (left) and nuclear (right) 
recoil distributions. 



measured a background- and statistics-limited ERC of 
7.6 X 10~^ between 52 and 110 keVr (1 contamination 
event). We observe no contamination events above 62 
keVr using the multibin method. 



IV. MAXIMUM LIKELIHOOD PSD 

We predict the PSD achievable in a large detector with 
many PMTs by assuming we can measure discrete times 
for each photoelectron detected in an event. We simu- 
late events in this detector using the PDFs measured in 
Sec. IIII Al and assume that the detection times of pho- 
toelectrons are measured without error. We develop a 
maximum likelihood PSD method and apply it to the 
simulated data. 

For an event that generates iV photoelectrons, we de- 
note the detection times as t = (ii, t2, ijv)- We define 
the log-likelihood function, InL, of this data following 
Ref. [Ill as 

N 

\nL = -m + J2Hng + Xf{t-to)), (17) 

i=l 

where 

TO = / {ng + Xf{t-toj)dt, (18) 

to is the time at which the energy deposit occurs, fit — to) 
is the PDF for the observed photoelectrons [Eq. |3] , A is 
the expected number of detected photoelectrons gener- 
ated by the event, T;, and Tg are the start and end time 
of the observation, and r^g is the background rate. 

In our simulations, we set rjjg — and to > T^. The 
maximum likelihood estimates of to and A are therefore 
ti, the first detection time, and A, where 



qF{T,-ti,Ti) + {l-q)F{T,-ti,T2) ' ' 

and F{T,t) = 1 — cxp(— T/t). Consequently, m = N. 
Additionally, the time window in our simulations is 40 ^s, 
so f{t — to)dt « 1. Thus, for our case the likelihood 
function of the observed data is well approximated as 

N 

i(t)«exp(-7V)iV^n (20) 

i=l 

Following Ref. [l^ and the discussion above, for each 
simulated event we determine L{t) for both the electronic 
and nuclear recoil event classes as modeled in Sec. IIII Al 
These values are denoted Le{t) and i„(t). We define the 
log-likelihood ratio. In R: 

lni? = lnie(t)-lni„(t). (21) 

Figure [m shows this log- likelihood statistic for simulated 
electronic and nuclear recoil events that yield 50 photo- 
electrons. The event is assigned to the nuclear recoil 
class if Ini? is less than an adjustable threshold that can 
be varied to increase the discrimination against electronic 
recoils at the cost of decreasing nuclear recoil acceptance. 
The Monte Carlo estimate of the distribution of In R has 
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prominent ripples even though the Monte Carlo estimates 
of L„ and L^, do not. In our Monte Carlo study, we sim- 
ulate events that yield a fixed number of photoelectrons 
according to PDFs with the assumption of perfect knowl- 
edge. Furthermore, we neglect dark current noise. In an 
actual experiment, we expect that imperfect energy res- 
olution due to variability in counting statistics and dark 
current effects would attenuate the ripples in the Ini? 
distribution. 




-"^O -30 -20 -10 10 20 30 40 
Log-likelihood statistic, In R 



FIG. 14: (Color online) Monte Carlo estimates of the log- 
likelihood statistics for events that yield 50 photoelectrons. 

To illustrate the maximum likelihood method, we sim- 
ulate a detector with a signal yield of 6 photoelec- 
trons/keVee, as might be possible in a detector with full 
PMT coverage. We neglect dark current and we set the 
discrimination threshold to accept 50% of nuclear recoils. 
For comparison, we simulate prompt fraction data by as- 
suming the number of prompt photoelectrons is a bino- 
mial random variable with an expected value determined 
by the measured PDFs, and the total number of photo- 
electrons is assumed to be known without error. This ide- 
alized binomial model predicts a much lower ERC than 
the statistical model discussed in Sec. IIIIBI 

The maximum likelihood discrimination method out- 
performs the idealized prompt fraction method, as shown 
in Fig. [ini We also expect that the maximum likelihood 
method will be more robust than the prompt fraction 
method to background noise. For the prompt fraction 
case, we select a threshold in prompt photoelectron space 
to yield a nuclear recoil acceptance probability as close 
to 0.5 as possible. Due to quantization effects, the ac- 
tual nuclear recoil acceptance probability varies about 
0.5, and there are sawtooth-like artifacts in the prompt 
fraction ERC curve shown in Fig. [13 

V. CONCLUSION 

Using a detector with a signal yield of 4.85 photoelec- 
trons/keVee, we measured the scintillation time depen- 



dence of electronic and nuclear recoils in liquid argon 




Energy (keVr) 

FIG. 15: (Color online) Predicted performance of maximum 
hkelihood and prompt ratio PSD methods for a detector yield- 
ing 6 photoelectrons/keVee (the energy axis has been scaled 
from keVee to keVr by use of a constant nuclear recoil scintil- 
lation efficiency of 0.29 as discussed in Sec. IIII Bl 



down to 5 keVee or 20 keVr. We developed a prompt frac- 
tion method of PSD in liquid argon, and for a nuclear re- 
coil acceptance level of 50%, we measured a background- 
and statistics-limited level of ERC to be 8.5 x 10~^ be- 
tween 52 and 110 keVr with no contamination events 
above 69 keVr. We also developed a multibin method 
of PSD, improving on the prompt fraction method by 
as much as an order of magnitude. With this method, 
we measured a background- and statistics-limited level 
of ERC of 7.6 X lO"'^ between 52 and 110 keVr for the 
same nuclear recoil acceptance of 50%. We modeled the 
observed prompt fraction data as the ratio of two nor- 
mally distributed, correlated random variables, where we 
assumed Np and Ni were uncorrelated; we discussed dis- 
crepancies between observed and predicted prompt frac- 
tion results. Finally, we developed a maximum likelihood 
method of PSD for a detector capable of measuring a dis- 
crete detection time for each observed photoelectron in 
an event. 
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